The FSSA algorithm now supports multivariate options allowing the user to perform multivariate functional singular spectrum analysis (mfssa) on multivariate functional time series (fts) objects. So, if we are tracking more than one family of curves that result from the same stochastic process, we can make our analysis richer by including these different types of information into the model. In order to illustrate the use of the FSSA algorithm on mutlivariate fts objects, we give a bivariate example with the Jambi data set. Satellite images of the Jambi region in Indonesia were taken every 16 days over the course of about 20 years using NASA’s MODerate-resolution Imaging Spectroradiometer (MODIS). The goal was to track changes in vegetation in the Jambi region by using pixel values (which are between 0 and 1) where a higher pixel value corresponds to more vegetation (greener areas) and a lower pixel value corresponds to less vegetation (less green areas). Two pixel indices were tracked known as the normalized difference vegetation index (NDVI) and the enhanced vegetation index (EVI) with globabl coverage at 250 meters squared. Our goal in this example is to use mfssa to uncover changes in NDVI and EVI data, which both came from the same process, in order to better detect changes vegetation in the Jambi region. This type of analysis can inform us of any patterns are present in the vegetation indices and if there is a general loss of vegetation in the region which would be reflected in the trend. Thinking about the analysis in a big picture, we expect to uncover seasonal changes in vegetation due to four season changes in a year. While this information is useful it would be more useful to detect any other strange changes in vegetation that might be explained by deforestation due natural or man-made causes.
Since density functions are more informative than using average or maximum pixel values, we start by converting each of our 448 images of pixel values into densities. In this sense, each density will give us the distribution of pixel values within that particular image taken in a day. This allows us to look at the probability density function of observing a pixel value between zero and one for both NDVI and EVI images for each of the 448 days that an image was taken. We estimate the density functions for both the NDVI and EVI data using the following.
require(fda)
require(Rfssa)
## Raw image data
NDVI=Jambi$NDVI
EVI=Jambi$EVI
time <- Jambi$Date
## Kernel density estimation of pixel intensity
D0_NDVI <- matrix(NA,nrow = 512, ncol = 448)
D0_EVI <- matrix(NA,nrow =512, ncol = 448)
for(i in 1:448){
D0_NDVI[,i] <- density(NDVI[,,i],from=0,to=1)$y
D0_EVI[,i] <- density(EVI[,,i],from=0,to=1)$y
}In this sense, we are forming functional data (fd) fd objects since the densities correspond to continuous random variables between zero and one. We have that one fd object contains functional elements that are the densities of the NDVI pixel values and the other contains densities of the EVI pixel values. We form these objects with the following code.
## Define functional objects
d <- 11
basis <- create.bspline.basis(c(0,1),d)
u <- seq(0,1,length.out = 512)
y_NDVI <- smooth.basis(u,as.matrix(D0_NDVI),basis)$fd
y_EVI <- smooth.basis(u,as.matrix(D0_EVI),basis)$fd
y=list(y_NDVI,y_EVI)At this point, we have two families of fd curves. One that corresponds to NDVI densities and the other that corresponds to EVI densities. These fd objects are really fts objects that were tracked over the course of 448 days. In this sense, we are really dealing with a multivariate fts object of pixel densities for both the NDVI densities and the EVI densities and so we convert the fd objects into a multivariate fts object using the following.
## Define functional time series
Y <- fts(y,time=time)We can visualize the multivariate fts object using line plots, heat maps, 3D surfaces, and 3D line plots. The line plots can be used to analyze patterns within a particular curve, the heat map can be used to find patterns between elements of the multivariate fts object, and the 3D plot options combine both types of information in one plot.
It is difficult to pick out any pattern within the curves of the mutlivariate fts object using the following line plot.
plot(Y,main="Density Estimations",xlab="Pixel Value",ylab=c("NDVI Pixel Density","EVI Pixel Density"), type="line")We get a little more information out of the heat map plot for the multivariate fts object as there might be some seasonality at play with regards to the four seasons in a year. This is to be expected as the probability of observing a green pixel in the colder months will drop and will also rise in the warmer months. We should note here that the first heat map corresponds to the NDVI images while the second corresponds to the EVI images.
plot(Y,main="Density Estimations NDVI and EVI",tlab="Time",xlab="Pixel Value", type="heatmap")The 3D surface and 3D line plots also show the seasons at play in influencing the probability of observing green pixels.
plot(Y, type = '3Dsurface', var=1,ylab = c("NDVI"),main = "Probability Kernel Density")plot(Y, type = '3Dline', var=1,ylab = c("NDVI"),main = "Probability Kernel Density")